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Abstract. We consider the problem of recovering an A-dimcnsional sparse vector x from 
its linear transformation y = Dx of M(< N) dimension. Minimizing the Zi-norm of a; under 
the constraint y = Dx is a standard approach for the recovery problem, and earlier studies 
report that the critical condition for typically successful ^-recovery is universal over a variety 
of randomly constructed matrices D. For examining the extent of the universality, we focus 
on the case in which D is provided by concatenating T = N/M matrices 0\, O2, ■ ■ ■ , Ot 
drawn uniformly according to the Haar measure on the M x M orthogonal matrices. By 
using the replica method in conjunction with the development of an integral formula for 
handling the random orthogonal matrices, we show that the concatenated matrices can 
result in better recovery performance than what the universality predicts when the density 
of non-zero signals is not uniform among the T matrix modules. The universal condition 
is reproduced for the special case of uniform non-zero signal densities. Extensive numerical 
experiments support the theoretical predictions. 



1. Introduction 

The recovery problem of sparse vectors from a linear underdetermined set of equations has 
recently attracted attention in various fields of science and technology due to its many 
applications, for example, in linear regression [I], communication [2], [3], [I], multimedia [5], 
[Ej, [7], and compressive sampling (CS) [8], [9]. In such a sparse representation problem, we 
have the following underdetermined set of linear equations 

y = Dx, (l) 

where y G M is an observation vector, D G M. MxN is the dictionary (or sensing matrix 
in CS), x G M. N is the sparse vector of interest, and M < N. Another way of writing (Tj[|) 
is that a large dimensional sparse vector x is coded/ compressed into a small dimensional 
vector y and the task will be to find the x from y with the full knowledge of D. For this 
problem, the optimum solution is the sparsest vector satisfying ([1]). Finding the sparsest 
vector is however NP-hard; thus, a variety of practical algorithms have been developed. 
Among the most prominent is the convex relaxation approach in which the objective is to 
find the minimum /i-norm solution to (TjQ). For the Zi-norm minimization, if x is i^-sparse, 
which indicates that the number of non-zero entries of x is at most K, the minimum K that 
satisfies ([1]) gives the limit up to which the signal can be compressed for a given dictionary 
D. An interesting question then arises: How does the choice of the dictionary D affect the 
typical compression ratio that can be achieved using the /i-recovery? 

Recent results in the parallel problem of CS, where D acts as a sensing matrix, reveal 
that the typical conditions for perfect /i-recovery are universal for all random sensing matrices 
that belong to the rotationally invariant matrix ensembles [10]. The standard setup, where 
the entries of the sensing matrix are independent standard Gaussian, is an example that 
belongs to this ensemble. It is also known that the conditions required for perfect recovery 
do not in general depend on the details of the marginal distribution related to the non-zero 
elements. On the other hand, we know that correlations in the sensing matrix can degrade 
the performance of Zi-recovery [H]. This suggests intuitively that using a sample matrix 
of the rotationally invariant ensembles as D is preferred in the recovery problem when we 
expect to encounter a variety of dense signals y. However, the set of matrix ensembles 
whose /i-recovery performance are known is still limited, and further investigation is needed 
to assess whether the choice of D is indeed so straightforward. 

The purpose of the present study is to fulfill this demand. Specifically, we examine 
the typical /i-recovery performance of the matrices constructed by concatenating several 
randomly chosen orthonormal bases. Such construction has attracted considerable attention 
due to ease of implementation and theoretical elegance [12], [13], [2] for designing sparsity 
inducing over-complete dictionaries for natural signals [15]. For a practical engineering 
scheme, audio coding (music source coding) [16] uses a dictionary formed by concatenating 
several modified discrete cosine transforms with different parameters. 
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By using the replica method in conjunction with the development of an integral 
formula for handling random orthogonal matrices, we show that the dictionary consisting of 
concatenated orthogonal matrices is also preferred in terms of the performance of Zi-recovery. 
More precisely, the matrices can result in better /i-recovery performance than that of the 
rotationally invariant matrices when the density of non-zero entries of x is not uniform 
among the orthogonal matrix modules, while the performance is the same between the two 
types of matrices for the uniform densities. This surprising result further promotes the use 
of the concatenated orthogonal matrices in practical applications. 

This paper is organized as follows. In the next section, we explain the problem setting 
that we investigated. In Section 3, which is the main part of this paper, we discuss the 
development of a methodology for evaluating the recovery performance of the concatenated 
orthogonal matrices on the basis of the replica method and an integral formula concerning the 
random orthogonal matrices. In Section 4, we explain the significance of the methodology 
through application to two distinctive examples, the validity of which is also justified by 
extensive numerical experiments. The final section is devoted to a summary. 



2. Problem Setting 



We assume that TV is a multiple number of M; namely, T = N/M = 2, 3, . . . G N. Suppose 
a situation in which an M x N dictionary matrix D is constructed by concatenating T 
module matrices Oi, 2 , . . . , Ot, which are drawn uniformly and independently from the 
Haar measure on M x M orthogonal matrices, as 

D= [O, 2 O t }. (2) 

Using this, we compress a sparse vector x° G M. N to y G M M following the manner of (IT]). 
We denote x° for the concatenation of T sub-vectors of M dimensions as 



x 



(3) 



\x T J 

yielding the expression 

y = Dx° = Oxxl + 2 x\ + ... + O T x T . (4) 

With full knowledge of O t (t = 1, 2, . . . , T) and y, the Zi-recovery is performed by solving 
the constrained minimization problem 

t 1 T 

^2 \\ x t\\i \ sub j- t0 V = ^20 t x t , (5) 



mm 



f=l 



t=l 
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where ||cc||i = X)t=i \ x i\ f° r ^ x = ^ anc ^ m i n x{/PO} generally denotes the 

minimization of f(X) with respect to X and G 1R (t = 1,2,... , T). At the minimum 
condition, {ai t } constitutes the recovered vector x in the manner of ([3]). 

For theoretically evaluating the Zi-recovery performance, we assume that the entries of 
x®, are distributed independently according to a block-dependent sparse distribution 

p t (x) = (1 - pt)5(x) + ptft(x), (6) 

where < pt < 1 means the density of the non-zero entries of the t-th block of the same size 
M in x° and ft(x) is a distribution whose second moment about the origin is finite, which 
is assumed as unity for simplicity. Intuitively, as the compression rate a = M/N = T" 1 
decreases, the overall density p = T" 1 J2t=i Pt U P t° which (jSJ) can successfully recover a 
typical sample of the original vector x° becomes smaller. However, precise performance may 
depend on the profile of p t . The above setting allows us to quantitatively examine how such 
block dependence of the non-zero density affects the critical relation between a and p for 
typically successful ^-recovery of x°. 



3. Statistical mechanics approach 

3.1. Statistical mechanical formulation 
Expressing the solution of (j5J) as 

x = lim / dxPp (x\y, {O t }) x, (7) 

/3^oo J 



where 



(x\y, {O t }) = Z~\fS- {x° t }, {O t })5 \y-J2 ex P ( -I 3 E H* 

= Z-\fi- {x° t }, {O t })5 (j2 Ot(x° t - x t )) exp [-P ll^lli 



,t=i / \ t=i 

and Z -1 (# {a#}, {O t }) = J l\J =1 dx t 5 fej=i ~ x tj) exp ^J =1 ||aj t ||iY consti- 

tutes the basis of our analysis. Equations ([7]) and (jHJ) mean that & can be identified with 
the average of the state variable x for the Gibbs-Boltzmann distribution (jSJ) in the vanishing 
temperature limit j3 — > oo. However, as (|HJ) depends on and {O t }, further averaging 
with respect to the generation of these external random variables is necessary for evaluating 
the typical properties of the Zi-recovery. Evaluation of such "double averages" can be carried 
out systematically using the replica method |17j . 

3.2. Integral formula for handling random orthogonal matrices 

In the replica method, we need to evaluate the average of Z n ((3; {x^}, {Ot}) for VnGN with 
respect to 0±, O2, ■ ■ ■ , Ot over the uniform distributions ofMxM orthogonal matrices. 
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However, this is rather laborious and is easy to yield notational confusions. For reducing 
such technical obstacles, we introduce a formula convenient for accomplishing this task 
before going into detailed manipulations. Similar formulae have been introduced for handling 
random eigenbases of symmetric matrices [T0l[T8l - l2T] and random left and right eigenbases 
of rectangular matrices [22|[23j. 



Let us assume that M-dimensional vectors u t 



x, 



1*1 



v 2 = M 1 1 rt 2 1 



characterized by their norms as V\ = M~ 
| it | denotes the standard Euclidean norm of the vector u. 
function 



x t (t 

..,V T 



- 1,2 
M _1 | 



, T) are 
2 , where 



For these vectors, we define the 



lim — In 

M->oo M 



lim — In 

M-s-oo M 



E°* w * 

t=i 



{Ot}, 



(9) 



/ (l£=i vo t 

where [f(X)] x generally denotes the average of f(X) with respect to X, and VO denotes 
the Haar measure of the M x M orthogonal matrices. Our claim is that by explicitly using 
v±, V2, ■ ■ ■ , Vt, © can be expressed as 

1 



extr 

Ai,A 2 ,...,At 



In 



,*=i 



t=i 



t=i 



t=i 



T 
2' 



(10) 



where extrx{/(^)} generally denotes the extremization of function f(X) with respect to 
X. Expression (jiop is derived from the fact that for fixed u t , O t u t moves uniformly on the 
surface of the M-dimensional hypersphere of radius \/Mv t when O t varies according to the 
uniform distribution of the orthogonal matrices; therefore, T>O t in (jH]) can be replaced with a 
spherical measure of M-dimensional vector du t 5(\u t \ 2 — Mv t ). For details, see Appendix A 



Function F({vt}) physically represents a characteristic exponent of the probability that 

= when 
2 = Mv t . 



M-dimensional vectors Ut (t — 1, 2, . . . , T) form a closed loop satisfying ^2 t=1 u t 



they are independently and isotropically sampled under the norm constraints of \u t 
For small T, the loop condition strongly restricts the region of {v t } to which F({v t }) is well 
defined. In concrete terms, F({vt}) diverges to minus infinity unless an equality 

vi = v 2 (11) 
and the triangular inequalities 

vi < v 2 + v 3 , v 2 < vi + v 3 , v 3 < vi + v 2 (12) 

are satisfied for T = 2 and 3, respectively. In particular, the constraint of (ITT]) requires us 
to deal with the system of T = 2 in a manner different from that of T > 3 except for the 
case of pi = p 2 , as discussed in the analysis of section 3.4. 
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3.3. Replica method 



Now, we are ready to apply the replica method for analyzing the typical property of the 
/i-recovery ([7j). For this, we evaluate the n-th moment of the partition function using the 
identity 

z n (P; R°}, {O t }) = / (f[f[dx a t 




K a=l t=l 

xexpf-^^^KHi), (13) 
\ =i t=i J 

which is valid for only tiGN. In the large system limit M — » 00, the rescaled logarithm of 
the moment, M _1 In [Z n {(3; {x^}, {Ot})] r x oy r 0t i, can be accurately evaluated for all n G N 
by using the saddle point method with respect to macroscopic variables Q® = M~ x \x^\ 2 , 
qf> = M~ l x^-x b t , and m" = M~ l x® -xf, where a,b — 1,2, ... ,n and t = 1, 2, . . . , T. Intrinsic 
permutation symmetry concerning the replica indices a — 1,2, ... ,n in ( JT3l) guarantees that 
there exists a saddle point of the form Q° = Q t , qf b = q t , and m" = m t , which is often 
termed the replica symmetric (RS) solution. As a simple and plausible candidate, we adopt 
this solution as the relevant saddle point for describing the typical property of the /i-recovery, 
the validity of which will be checked in section 3.5. The detailed computation is carried out 
as follows. 



3.3.1. Energetic part Let us consider averaging ( fl3l) with respect to {Ot} and define for 
each fixed set of {x^}: 

11 / T 

II' 5> 



'tu t 



a=l \t=l 



(14) 



{Ot} 



where u" = x® — x%. When {x®} is placed in the configuration of the RS solution, the 
expression 

/ MR t Mr t ■ ■ ■ Mr t \ 



[u\ u\ ... <] T x [u] u\ ... <] = 

/ M(R t -r t + nr t ) 



Mr t MR t 
\ Mr t Mr t 



E x 







M(R t - r t ) 



MR t J 

\ 


" x E 1 



V 



(15) 



• • • M(R t - r t ) J 

holds for each t, where T stands for matrix transpose, Rt — Qt — 2m t + pt, and r t = 
qt — 2m t + pt- E = [e\ e-i ■ ■ ■ e„] denotes annxn orthogonal matrix composed of the vector 
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ex = (n l ' 2 , n l / 2 , . . . , n 1 / 2 ) T and an orthonormal set of n — 1 vectors e 2 , e 3 , . . . , e n that 
are orthogonal to e±. This indicates that [u\ u 2 . . . it™] may be expressed as 

[uj u 2 . . . <] = [fij u 2 . . . <] x £ T , (16) 

by using a set of n orthogonal vectors {ill }> whose norms are given as \u]\ 2 = M(R t — r t + 
nr t ) = M(Qi-g 4 + n(gt-2m t + p t )) and \ii a t \ 2 = M{R t -r t ) = M(Q t -q t ) for a = 2, 3, . . . , n, 
along with an n x n orthogonal matrix E that does not depend on t. This guarantees the 
equality X({w^}) = I({u^}). Furthermore, condition n«M and the orthogonality of {ii^} 
among the replica indexes a — 1,2, ... , n allows us to evaluate the average concerning {O t } 
independently for each index a when computing I({ii^}). This, in conjunction with 0, 
provides each set {x®} of the RS configuration with an expression of (fT4|) as 

± lnX(K» = F({Q t -q t + n(q t - 2m t + p t )}) + (n - l)F({Q t - q t }). (17) 

The right hand side of (TTTT) is likely to hold for n e M as well, although (fT4|) is defined 
originally for only fieN. 



3.3.2. Entropic part On the other hand, inserting identities 1 = J MdQtS(\x1 1 2 — MQt), 
1 = / Mdq t 8{x1 ■ x\ - Mq t ) and 1 = f Mdm t S(x° t ■ x a t - Mm t ) (a,b = 1, 2, . . . , n; t = 
1,2, ... ,T) into (fT3|) and taking an average concerning {a^ }, in conjunction with integration 
with respect to dynamical variables {x^ }, result in the expression 



V({Q t ,qt,m t }) 



T n 



T n 



niT* 

,t=l a =l 



1 exp 



T 



\5(xt ■ x\ - Mq t 



x 



,t=l a>b 



T 

n 

t=i 



Hl[5(\x*\ 2 - MQ 

t=l a =l 

~[S(x° t ■ x a t - Mm t 



a=l 




dQexp [M{A{{Q t ,q t ,m t },{Q a t ,q a t b ,m a t }) + B{{Q a t ,qf ,m a t }) 

for a fixed set of {Qt, qt, mt}- The conjugate variable Qt is introduced for expressing a 
delta function as 5{\x a t \ 2 - MQ t ) = (An^/^iy 1 f+^™dQ%exp (-(Q?/2) ( | ccf | 2 - MQ t )), 
and similarly for qf b and m". Notation dQ stands for an integral measure 
lL=i ( lla=i dQt Yl a >b dqt b lla=i dm? ) ) an d the functions on the right hand side are defined 



as 



T 



A({Qu ft, m t }, {4°, gf , m?}) = £ ( £ 9l9* _ £ gf ft _ £ ™ >t ) (19) 



and 



5({4,gf,<}) = ^ln 



i=l 



t=l \a=l 



\dxl 



\a=l 



a>b 



o=l 



exp 



Efiab a b 
qt ^t^t 



a>b 



a=l a=l / J . 

where the average for x\ is taken according to ([6]). The expression of (fl8l) indicates that its 
rescaled logarithm is accurately evaluated using the saddle point method with respect to the 
conjugate variables in the large system limit M — > oo. In addition, the replica symmetry 
guarantees that the relevant saddle point is of the RS form as Q1 = Q t , q1 b = qt, and 



m" = rh t . As a consequence, the evaluation yields 



— \nV({Q t ,qt,m t }) 



+ ln 



Dz\ I dxexp 




extr 

t,qt,m t 



nQtQt 



n(n 



l)<Mt 



nm t m t 



(Qt + QtW 



+ WGt z + ™,t x t ) x ~~ ft I x I 



,(21) 



where Dz = dzexp(—z 2 /2)/\/2ir denotes the Gaussian measure. This is also likely to hold 
for net, although (1181) is originally defined for only n G N. 



3.3.3. Free energy and saddle point equations The replica method uses the identity 



■(/3M)" 1 [lnZ(/3;K},{O t })] 



{*?},{Ot} 



lim n ^ (9/9n)(/3M)- 1 In {[Z n {fo {x° t }, {O t } 



)]| a .o} {o t }j f° r evaluating the typical free energy density. The above argument indicates 

that M -1 In ([Z n (f3; {x®}, {Ot})]{ x oy {o t }) can De computed by extremizing the sum of ffTTj) 
and (12TI) with respect to {Q 4 ,g 4 ,m t }. Furthermore, the obtained expression is likely to hold 
for n 6 M, although the calculations are based on ( Fl3|) that is valid for only n G N. We, 
therefore, take the limit of n — > utilizing the expressions of (IT7|) and ( 121]) for n G [0, 1] as 
well. In particular, in the limit of (3 — > oo, which is relevant in the current problem, the 
expression of the free energy of the vanishing temperature is expressed as 

— [^m{* t h{o t })] {x o h{0t} 

f dF({ Xk }) 

dxt 



— lim 

T 



t=i i 



(Q t -2m t + pt) + 



QtQt 



XtXt 



+ rhtmt 



Dz 



where 



( \/xtZ + m t x°;Qt 




X t z + m t x I x + \x 



(22) 



(23) 



and rescaled variables are introduced as /3(Q t — qt) — > Xt, (jt/P 2 Xt, (Qt + Qt)/P — > Qt, and 
rh t / (3 tth to properly describe the relevant solution in the limit of /3 — > oo. Extremization 
is to be performed with respect to {Q t , Xu m t, Qt, Xt, rht}- 
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Similar to earlier studies, at the extremum characterized by a set of the saddle point 
equations 

,dF({ Xk }) = fit 

Xt 



Qt 



rh t = -2 



dxt 
d 2 F({ Xk }) 
\ dxtdxs 

dF{{ Xk }) = Rt 
Xt 



(Q a - 2m s + p a ), 



Qi 



Xt 



m t 



Dz 
Dz 

Dz 



dxt 

X 2 (^/x~tz + m t x°; Q t ) 

dX 2 (Vx~tZ + m t x° t ;Q t ) 
d(Vx~tz) 

%t x (Vx~tz + rn t x° t ]Q t ) 



(24) 
(25) 
(26) 

(27) 
(28) 
(29) 



Q t and m t (t = 1,2, ... ,T) physically denote the macroscopic averages of the recovered 
vector x t as M~ l \\x t \ 2 }^ y | 0fe j. and M" 1 [x° t ■ x t ]y x oj y Qk y, respectively Here, R t = 

A^ 1 / (j2f = i^t lS ^ is provided by the extremum solution of ([IP]) for {v t } = {xt}, and 

r (h-i)/Q, h>i, 

X{h;Q) = -d<p(h;Q)/dh = I 0, \h\ < 1, (30) 

[ (h + l)/Q, h<-l. 

For T = 2, a Lagrange multiplier should be exceptionally introduced in (125)) for enforcing 
Xi = X2- 



3.4- Critical condition for l\ -recovery 

The success of the Zi-recovery is characterized by the condition in which Q t = m t = pt 
is satisfied at the extremum for Vt = 1,2, ... ,T. Therefore, one can evaluate the critical 
relation between a = 1/T and p = T^ 1 Ylt=i Pt by examining the thermodynamic stability of 
the success solution Q t = m t = pt (t = 1, 2, . . . , T) for the saddle point equations (|24l)-(f29l) . 

We assume T > 3 for a while, since an exceptional treatment is required for T = 2. 
For obtaining the success solution, it is necessary that Xt = and Q t = rh t = +oo hold. 
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Expanding ( )27|) -( l29|) under the assumption of Q t = rh t ^> 1 yields 

Q t ~ Pt - -^j- + (2(1 - Pt )G(xt) + Pt(xt + 1)) + 0(4" 3 ), (31) 



27iQ t Q 

m t ^p t -J^ + 0{Qf), (32) 
V lixQt 



xt H 2 ' 1 - ft » Q fe" /2 » + '", (33) 



where Q(x) = Dz 



, p -l/(2x) 

G(x) = (x + l)Q(x~ 1/2 ) -x 1/2 , (34) 

V 27T 



and we used [(x°) 2 ] x o = p i; which is derived from the assumption (jBJ). These result in the 
expression 

n o i 2(1 -p f )G(xf) + Pt(Xf + 1) , n /A-3N 

Q t -2m t + p t ~ \-0{Q t ). (35) 



Q 2 



In addition, (fTUj) indicates that 



<9 2 F({ Xfc }) 1«9A S , 1 1 /c^V* , 1 



a***. 2 ^ + 2^ ~ 2 as: ) ts + M 6ts m 

holds, where {A t } is the solution of T coupled equations 

"-hi 1 - (37) 

Differentiating this yields 

dXm (1 — 2R m ) R m Rn 



(Jmn + K mn ), (38) 

where we denote J = (J mn ) = ((1 - 2R m )K^5 mn ) and K = (K mn ) = ((R m R n ) (A m A„) _1 ) . 
This expression makes it possible to evaluate (dA m /dx n ) = (^Xm/^A n ) _1 as 



= -^2(-l) p (J- 1 K) p J- 1 = -T(I + M)- 1 T, (39) 

p=0 

where / = (S mn ), T = (A m (l - 2R m )- 1 / 2 5 mn ) and M = {R m R n {l - 2R m )- l ' 2 {l - 2i? n )" 1 / 2 ). 
Inserting (1351) . (136|) . and (|39|) into (f25|) yields a set of equations to determine {xt} for a given 
set of non-zero densities {pt} 

2{l-p t )G{xt)_+Pt{Xt + l) 
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- J2(S(I + M)~ 1 S) ts (l - R s ) 2 (2(1 - p s )G(x s ) + Ps(Xs + 1)) , (40) 

s=l 

where we set S = (S mn /(R m y/l — 2R m )) and used the relation T/Q t = (1 — Rt)S, which is 
obtained from (124|) and (1371). 

Equations (125]) and ( 133|) indicate that the critical condition for making Xt — stable is 
expressed as 

R t = 2(l-p t )Q(x; 1/2 )+Pt. (41) 
Furthermore, the condition 

f> = l (42) 
t=i 

must hold by the definition of R t . 

For characterizing the critical condition for the success of the /i-recovery, let us suppose 
that the set of non-zero densities {pt} is provided as a function of a single parameter p as 
{pt(p)}. For T > 3, the critical situation is specified by an appropriate set of 2T + 1 variables 
of p, {xt}, and {Rt}. These are provided by 2T + 1 conditions of fT40|) -f l42|) . 

On the other hand, the critical condition for T = 2 is provided differently from 
( )40|) -( 142|) because the constraint of Xi — X2 for keeping F({xi,X2}) well defined requires 
Ri = R2 = 1/2 for any pair of p\ and p2- Explicitly, the condition is provided by the 
following four coupled equations: 

Xi = 4(1 - Pi(/i))G(xi) + 4 Pl (/i)(xi + 1) + V, (43) 
X2 = 4(1 - p 2 (^))G(x 2 ) + 4p 2 (p)(x 2 + !)-??, (44) 
2(1 - Pl (/i))Q(xr 1/2 ) + 2pi(/i) = \, (45) 

2(1 - p 2 (p))Q(x 2 " V2 ) + 2p 2 (/i) = \, (46) 

where 77 is a Lagrange parameter for enforcing xi = X2- These determine four variables of 
Xi; X2, f]-, P at the critical condition. 

Equations p3]) -p2" j) for T > 3 and (i3])-(jl6]) for T = 2 constitute the main result of 
this paper. For T = 2, an equivalent result can also be obtained in a slightly different 
manner 241. 



3. 5. Validity of RS evaluation 

The above calculation can be generalized to arbitrary levels of replica symmetry breaking 
(RSB). For example, under one-step RSB (1RSB) ansatz, where n replicas of each block t 
are classified into n/x groups of an identical size x and their overlaps are assumed to be 
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xl ■ x b t /M = q lt if a and b belong to the same group and got(< q\t) otherwise, (fT7|l and (}2~Tj) 
are modified as 

lnl({<}) = F({Q t - q lt + x(q lt - q ot ) + n(q ot - 2m t + 
+ g - l) F({Q t - g u + x(g lt - g ot )}) + (n - ^) F({Q t - q lt }) (47) 



and 



T 



1, ,,/r n u , ] nQ t Q t n(x -l)(q lt q lt -q ot q ot ) 

— tnV{{Qt,qit,qot,rn t }) = . extr 



^ i=1 Qt,qit,qat,rh t 



n(n - l)q ot q t . , , 
nm t m t + In 



/Dzi /cfeexp(<l>) V , (4 




respectively, where $ = -(Qt + qit)x 2 /2 + (Vgit = QotZi + v/got^o + m t xf)x - (3\x\. 

In the 1RSB framework, the RS solution is regarded as a special solution of At = 
qit — qot = and At = g\ t — qot = (t = 1, 2, . . . , T). In the limit of n — > and /3 — > oo 
keeping x = /3(Qt — git) ~ 0(1), the critical condition that a solution of A t > and 
A t > bifurcates from the RS solution, which corresponds to the de Almeida-Thouless (AT) 
condition [23] of the current system, is expressed as 

det (i - HH^j = 0, (49) 

where H = (2(d 2 /d X td X s)F({Xk})) and H = (5 ts j Dz[(dX(Vxt^+mx° t ; Q t )/d(Vx~tz)) 2] A 
For T > 3, 



j j. 



rrrr i ( /+M)r,'(i -fl,y , (50) 



holds for the success solution at the critical condition of the Zi-recovery. Equation ( 1501) yields 
an eigenvalue of unity whose eigenvector is given as t>i oc ((1 — i? t ) _1 ), which makes fT4*9|) 
hold. Similarly, ( 149]) is also satisfied at the critical condition of the /i-recovery of T = 2. 
These validate our RS evaluation in terms of the local stability analysis, although further 
justification with other schemes, such as comparison with numerical experiments, is necessary 
for examining possibilities that the RS solution becomes thermodynamically irrelevant due 
to discontinuous phase transitions. 

4. Case studies 

Let us examine the significance of the developed methodology by applying it to two 
representative examples. 
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Table 1. Comparison of critical values of total non-zero density p = T^ 1 ^ Jt=l pt for 
uniform and localized densities. The values are rounded off to the fourth decimal. The 
values for uniform density is identical to those of the rotationally invariant matrix ensembles 
for compression rate a = M/N = 1/T. 



T 



2 3 4 5 6 7 8 



uniform 
localized 



0.1928 0.1021 0.0668 0.0487 0.0378 0.0308 0.0257 
0.2267 0.1190 0.0780 0.0566 0.0438 0.0354 0.0294 



4-1. Uniform and localized densities 

We consider the uniform density case of pt — p {t — 1, 2, . . . , T) as the first example, where 
the uniformity allows us to solve fl24^) -f )29|) setting all variables to be independent of t as 
Xt = X- I n particular, setting A t = A simplifies the expressions of fl24^) -f )26|) providing 
R t = 1/T and 2(<9 2 / dxtdXs) F ({xk}) = X~ 2 - This makes it unnecessary to deal with the 
saddle point problems of T = 2 in an exceptional manner. As a consequence, the critical 
condition of the /i-recovery is expressed compactly by using a pair of equations as 



for both T > 3 and T = 2. By setting a = 1/T, these provide a critical condition identical to 
that obtained for the rotationally invariance matrix ensembles in earlier studies [T0 | |2"6" | |2"7] . 
This indicates that for vectors of the uniform non-zero density, the /i-recovery performance 
of the concatenated orthogonal matrices is identical to that of the standard setup provided 
by the matrix of independent standard Gaussian entries. 

However, this is not the case when the non-zero density is not uniform. As a distinctive 
example, we examined the case of localized density, which is characterized by setting pi = Tp 
and pt = for t = 2, 3, . . . , T. Table [Hand Figure [1] show critical values of the total non-zero 
density p = T" 1 Ylt=i Pt given the compression rate a = 1/T for the uniform and localized 
density cases. These show that the concatenated matrices always result in better Zi-recovery 
performance for vectors of the localized densities, and the significance increases as T becomes 
smaller while matrices of rotationally invariant ensembles result in identical performance as 
long as p is unchanged. This indicates that, in addition to their ease of implementation and 
theoretical elegance, the concatenated orthogonal matrices are preferred for practical use in 
terms of their high recovery performance for vectors of non-uniform non-zero densities. 

4-2. Numerical justification 

To justify our theoretical results, we conducted extensive numerical experiments of the l\- 
reconstruction. Figures [2] (a) and (b) depict the experimental assessment of the critical 



X 



T(2(l-p)G(x)+p(x+l)), 



(51) 
(52) 
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Figure 1. (Color online) Critical values p c of total non-zero density versus compression 
rate a. Circles and crosses correspond to the localized and uniform densities, respectively, 
for a = 1/2, 1/3, ...,1/8. The curve represents the relation between p c and a for the 
rotationally invariance matrix ensembles. Crosses coincide with values of the curve for 
a = l/T (T = 2,3,...). 

threshold for T = 2 and T = 5, respectively. The case of an i.i.d. standard Gaussian 
dictionary is also plotted for comparison. Given fixed values of T and N, a trial was started 
with an empty vector x° and a concatenated orthogonal dictionary generated from a set 
of T standard i.i.d. M x M Gaussian matrices using QR-decomposition. Based on the 
relative densities {pt}, one sub-vector x® was then randomly chosen and assigned a non- 
zero component drawn from the standard Gaussian ensemble. Matlab algorithm "linprog" 
from Optimization Toolbox was used to solve the Zi-minimization problem and obtain the 
reconstruction x. The reconstruction was deemed to be a success if ||cc° — x\\x < 10 -6 and a 
failure otherwise. Given a successful reconstruction, we again randomly chose one sub- vector 
x® based on the densities {pt} and inserted a non-zero component drawn independently from 
the standard Gaussian ensemble into it. The process was continued until the original vectors 
{x®} had {K t } non-zero components and the reconstruction was deemed a failure, that is, 
\\x° — x\\i > 10~ 6 . The critical value K c = ^2 t K t — 1 was recorded and the experiment 
was started again using a new independent dictionary and an empty vector x°. For each 
value of T and N, we carried out 10 6 independent trials. The experimental critical density 
was defined as p c (T,N) = K c /N, where 7TT denotes the arithmetic average over the trials. 
For all system sizes, we also computed the experimental per-block densities {K t /M} and 
checked that they were close to the desired densities {pt} after the 10 6 trials. For fixed 
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0.005 0.01 
1 / N 



Figure 2. (Color online) Experimental assessment of critical densities for ^-reconstruction. 
Experimental data (see the main text) were fitted with a quadratic function of 1/N and 
plotted with solid lines. Concatenated orthogonal "O" and i.i.d. standard Gaussian 
"G" basis under uniform and localized densities. Filled markers represent the predictions 
obtained through the replica analysis. Extrapolation for N = TM — > oo provides the 
estimates for the critical values p c . (a) T = 2, where the markers correspond to simulated 
values M = 8, 9, . . . , 25, and (b) T = 5, where the markers correspond to simulated values 
M = 13, 14, ... , 30. 



Table 2. Comparison between experimental and theoretical assessments of critical values 
of total non-zero density p c for uniform and localized densities. The values are rounded off 
to the fourth decimal. For all cases in the table, the differences between the two estimates 
are most evident in the last digits. The excellent agreement to the experimental assessments 
justifies our theoretical treatment. 



T 


2 


3 


4 


5 


uniform (experiment) 
uniform (theory) 


0.1927 
0.1928 


0.1019 
0.1021 


0.0670 
0.0668 


0.0487 
0.0487 


localized (experiment) 
localized (theory) 


0.2264 
0.2267 


0.1196 
0.1190 


0.0779 
0.0780 


0.0567 
0.0566 



T, the experimental data points p c (T,N) were fitted with a quadratic function of 1/N. 
Extrapolation for N — > oo provided the experimental estimates of the critical densities, as 
listed in Table |2] in which the theoretical estimates in Table [T] are also listed for comparison. 

Comparing the theoretical and experimental results confirms the accuracy of the replica 
analysis. From Figure [2} we observe that for the finite-sized systems, the T orthogonal 
dictionaries seem to always provide higher thresholds p c than the Gaussian one, even for 
uniform densities. 



15 



5. Summary 

In summary, we investigated the performance of recovering a sparse vector x G M. N from 
a linear underdetermined equation y = Dx £ IR M when .D is provided as a concatenation 
of T = N/M independent samples of M x M random orthogonal matrices and the l\- 
recovery scheme is used. Performance was measured using a threshold value of the density 
p of non-zero entries in the original vector a? , below which the /i-recovery is typically 
successful for given compression rate a = M/N = T~ l . For evaluating this, we used the 
replica method in conjunction with the development of an integral formula for handling the 
random orthogonal matrices. Our analysis indicated that the threshold is identical to that 
of the standard setup for which matrix entries are sampled independently from identical 
Gaussian distribution when the non-zero entries in x° are distributed uniformly among T 
blocks of the concatenation. However, it was also shown that the concatenated orthogonal 
matrices generally provide higher threshold values than the standard setup when the non-zero 
entries are localized in a certain block. Results of extensive numerical experiments exhibited 
excellent agreement with the theoretical predictions. These mean that, in addition to their 
ease of implementation and theoretical elegance, the concatenated orthogonal matrices are 
preferred for practical use in terms of their high recovery performance for vectors of non- 
uniform non-zero densities. 

Promising future studies include performance evaluation in the case of noisy situations 
and development of approximate recovery algorithms suitable for the concatenated 
orthogonal matrices. 
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Appendix A. Derivation of (1101) 

When O t are sampled independently and uniformly from the Haar measure of the M x M 
orthogonal matrices, Ut = Otu t (t = 1,2, ... ,T) are distributed independently and uniformly 
on the surfaces of the M-dimensional hyperspheres of radius \fMv t for a fixed set of M- 
dimensional vectors u t satisfying |n t | 2 = Mvt- This means that the integral of (EJ) can be 
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expressed as 

I (l£i VO t ) 5 (Ef =1 O t u t 



/ (nf=i du t 8(\u t \ 2 - Mv t )) 5 (Ef =1 ut 
l{Ut=i du t 5(\u t \ 2 - Mv t 



(A.l) 



We insert the Fourier expressions of <5-function 

-loo 



i /"+V-100 /A \ 

5(\u t \ 2 - Mv t ) = j^j= J ^_ dA t exp (^-y (N 2 - Mv t )j (A.2) 

and 

' (t *) - <^7 £ * -p (^ft • g «,) ) (A.3) 

into the numerator of (1A.1|) . and carry out the integration with respect to {ut}, where 
k = (k\, /c 2 , . . . , This yields the expression 



oc 



oc 



,t=i / \t=i 



U^7 M ' 2 expfe^ fe^T 1 • 



2 

,t=i / \t=i / \t=i 



Evaluating this by means of the saddle point method with respect to A t (t = 1,2, ... ,T) 
results in 



lim — 



h la (./ fn^iw 2 - Mv <n s (fj"<)) 

\t=i / t=i t=i 



Similarly, the denominator of ( 1A.1I) is evaluated as 



lim — In 

M->oo M 



Substituting (tOjl . (TA~5l) . and (fA~6|) into US) leads to the expression of ( ITUl) . 
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